
% Function that computes standard errors with block bootstrap

function [se_theta,se_r1,se_r2,se_coef_ri_grp1,se_coef_rc_grp1,se_coef_ri_grp2,...
    se_coef_rc_grp2] = stderr(data,scf,k,l,m,sa1,invV,mhat_scf_param4)

    % Inputs
    % data      : data matrix
    % scf       : SCF data matrix
    % k         : number of model parameters: 4
    % l         : block length, i.e., number of years in a block
    % m         : number of bootstrap iterations
    % sa1       : fixed income asset share of group 1
    
    % Extract log rate variables
    ln_ri       = data(:,4);
    ln_rc       = data(:,5);
    
    % Number of blocks
    nb          = size(data,1)-l+1;
    
    % Construct (overlapping) bootstrap blocks
    blk         = cell(nb,1);
    for i = 1 : nb
        blk{i}      = data(i:(i+l-1),:);
    end

    % Pre-allocate matrices for bootstrap estimates
    theta       = zeros(k,m);
    r1          = zeros(size(data,1),m);
    r2          = zeros(size(data,1),m);

    % Loop over bootstrap iterations
    for j = 1 : m
        
        % Display bootstrap iteration
        fprintf('Bootstrap Iteration #%d\n',j);

        % Extract sample with replacement
        rng(j+1);
        blkboot     = datasample(blk,nb);
        bsample     = vertcat(blkboot{:});

        % Extract bootstrap sample variables
        y1          = bsample(:,1);
        y2          = bsample(:,2);
        a           = bsample(:,3);
        ri          = bsample(:,4);
        rc          = bsample(:,5);
        
        % Vector of sample moments
        moms        = mhat(y1,y2,a,ri,rc);
        
        % Compute sample moments from SCF in the 12-parameter model
        if k==4
            moms_scf    = mhat_scf_param4(9:14,1);
        end
        
        % Estimate parameter vector
        % Compute model-implied, group-specific rates of return
        if k==4
            theta(:,j)          = cmd_ipopt4(moms,moms_scf,sa1,invV,mhat_scf_param4);
            [r1(:,j),r2(:,j),coef_ri_grp1(:,j), coef_rc_grp1(:,j), coef_ri_grp2(:,j), coef_rc_grp2(:,j)] = computerates(theta(:,j),4,ln_ri,ln_rc);    
        end
        
    end
    
    fprintf('dimension of coefficient bootstrap')
    size(coef_ri_grp1)
    
    % Compute standard deviations of bootstrap estimates
    se_coef_ri_grp1 = std(coef_ri_grp1,[],2);
    se_coef_rc_grp1 = std(coef_rc_grp1,[],2);
    se_coef_ri_grp2 = std(coef_ri_grp2,[],2);
    se_coef_rc_grp2 = std(coef_rc_grp2,[],2);
    
    se_theta    = std(theta,[],2);
    se_r1       = std(r1,[],2);
    se_r2       = std(r2,[],2);
